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FINITE  DIFFERENCE  CALCULATIONS  FOR  HYDRODYNAMIC 
FLOWS  CONTAINING  DISCONTINUITIES 

ABSTRACT 

In  this  paper  we  show  how  to  calculate  the  steady 
hypersonic  inviscid  flow,  including  a  detached  shock, 
around  a  blunt  body.   The  steady  flow  is  obtained  as  the 
limit  for  large  time  of  time  dependent  flow,  starting 
with  plane  flow  impinging  on  the  body.   The  transient 
flow  is  the  solution  of  a  mixed  initial-boundary  value 
problem  for  the  partial  differential  equations  of 
inviscid  fluids  which  is  solved  by  a  difference  scheme 
proposed  by  Lax  and  Wendroff.   Our  calculations  show 
that  by  itself  this  difference  scheme  tends  to  be 
unstable  and  does  not  converge  to  the  steady  flow; 
by  adding  an  artificial  viscosity  term  we  have  succeeded 
in  stabilizing  the  calculation.   In  Section  4  we  give  a 
fairly  convincing  theoretical  explanation  of  this 
stabilizing  effect  and  derive  a  new  stability  condition. 

Both  plane  and  cylindrical  symmetries  are  considered; 
in  the  cylindrical  case  we  use  a  variant  of  Richtmyer's 
twostep  version  of  the  Lax-Wendroff  difference  scheme. 
This  method,  as  does  Richtmyer's,  requires  much  fewer 
arithmetic  operations  as  compared  to  the  one-step  method. 
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NYO-1480-33 

FINITE  DIFFERENCE  CALCULATIONS  FOR  HYDRODYNAMIC 
FLOWS  CONTAINING  DISCONTINUITIES 

Samuel  Z.  Burstein 

1 .   Introduction 

This  report  continues  the  exploration  of  certain 
explicit  difference  methods  for  fluid  dynamic  computations 
In  three  Independent  variables.   The  application  of  the 
Lax-Wendroff  method  was  extended  to  a  more  difficult 
problem;  the  detached  shock  problem.   The  method  of 
solution  Is  to  consider  the  problem  to  be  of  the  mixed 
type,  an  Initial  and  boundary  value  problem.   By 
prescribing  boundary  conditions  both  on  the  body  and 
at  the  upstream  and  downstream  regions  and  by  prescribing 
the  Initial  data,  a  time-dependent  solution  was  generated 
from  the  above  conditions.   Solutions  for  large  times  were 
obtained  and  these  correspond  to  the  steady  state.   Both 
plane  and  axlsymmetrlc  flows  were  considered  and  several 
difference  methods  were  tested  and  compared. 

In  a  previous  paper  [1]  results  were  presented  for 
a  series  of  calculations  relating  to  compressible  flows 
In  channels.   These  flows  contain  oblique  shocksand  for 
some  shapes  Mach  reflections.   It  was  noted  that  If  Mach 
reflection  took  place,  there  were  oscillations  present 
In  the  downstream  region  of  the  flow  field.  I.e.,  behind 
the  normal  shock  wave  and  extending  to  the  boundary;  no 
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such  oscillations  were  detected  for  the  case  of  regular 
reflection.   It  might  be  thought  that  the  difference 
between  the  two  cases  can  be  tentatively  explained  as 
follows:   the  treatment  of  the  downstream  boundary  by 
extrapolation  introduced  errors.   In  subsonic  flow  these 
errors  propagate  upstream  and  affect  the  accuracy  of  the 
solution  but  are  swept  out  of  the  mesh  in  supersonic  flows. 

However  severe  difficulties  were  encountered  in  the 
numerical  calculation  of  the  detached  shock  problem.   The 
solution  of  these  difficulties  shed  light  on  the  behavior 
of  the  difference  equations  and,  in  return,  on  the  nature 
of  the  oscillations  encountered  in  the  Mach  reflection 
calculation . 

The  author  would  like  to  acknowledge  that  during  the 
course  of  this  investigation  he  had  many  helpful  discussions 
with  Professor  Peter  Lax  --  all  of  which  made  this  work 
possible.   I  also  thank  Professor  Robert  Richtmyer  for  his 
comments  on  the  paper. 

2.   Hydrodynamic  differential  equations 

The  conservation  of  mass,  momentum  and  energy  in 
inviscid  time-dependent  fluid  motion  are  expressed  by 
equations  of  divergence  form 

(2.1)  div  (f)  =  0  . 
Here 

(2.2)  ^  =  (w  ,  f(w),  g(w)) 

and  div  is  the  partial  operator  in  three  dimensional 

space-time 
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(2.3) 
Hence 


dlv 


n,^  +  (  ),^  +  (  ), 


y 


^'t  +  ^'x  + 


y 


0  ; 


w  = 


w,  f  and  g  are  vectors.   The  components  of  w  are  mass, 
momentum  In  the  x-  and  y-directlon  and  the  total  energy 
all  per  unit  volume 

pu 
pv 

The  vectors  f  and  g  are  non-linear  vector  valued  functions 
corresponding  to  the  fluxes  in  the  x-  and  y-directions  of 
the  quantities  w. 

For  the  case  of  axially  syiranetric  flows  the  equations 
of  motion  take  the  slightly  different  form. 


(2.4) 


div  (j)  =  h 


V; 


T 


where  the  vector  h  =  -  — (w  +  tt),  tv     =    (0,0, 0,p)  and  the 
coordinate  y  stands  for  the  radial  distance  from  the  axis 

of  symmetry.   Equation  (2.4)  may  be  written  in  almost 

T 
conservation  form  if  the  vector  Jl      =  (0,0, p,0)  is 

introduced  by 

0  \ 


(2.5) 


-  yh  + 


0 


=  g 


Voy 


Expanding  Eq .  (2.4)  by  using  the  relations  (2.2),  (2. 5), 
(2.4)  and  (2.5)  yields 
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C2.6) 


w,,  +  f,^  +g,y  =  -  i  . 


Thus  the  vector  ^   is  replaced  by  y^   in  Eq .  (2.4)  while  the 
inhomogeneous  term  £    contains  only  one  nonzero  component 
in  the  radial  momentum  equation  compared  to  four  contained 
in  h.   For  net  points  not  on  the  axis  of  symmetry,  Eq .  (2.6) 
should  be  used  for  the  generation  of  the  difference  scheme 
in  axially  symmetric  flows. 

The  pressure  p  may  be  expressed  in  terms  of  the 
conservation  variables  and  the  ratio  of  specific  heats  of 
the  gas,  7,  through 


(2.7) 


(7-  1)  [E-(m2  +  n2)/2p) 


where  m  =  pu  and  n  =■  pv.   With  the  use  of  Eq.  (2.7),  the 
vectors  f(w)  and  g(w)  become 


f(w)  = 


g(w).  - 


m 


mn 

P 

(1-7)    /■  2  ,   2,  ,    mE 
^  2   m  (m  +  n  )  +  7  — 


n 

nm 
P 


For  axlally  symmetric  flows  p  Is  replaced  by  yp 
everywhere  in  the  above  expressions. 

5.   Difference  schemes  of  second  order  accuracy 

The  difference  scheme  that  Lax  and  Wendroff  used 
for  one-dimensional  calculations  [3]  and  the  two 
dimensional  version  which  was  analyzed  in  [2]  is  based 
on  the  Taylor  expansion 

C3.I)   w(t  +  At)  =  w(t)  +  At  w,^  +  At^/2  w,^^  . 

The  error  in  the  above  expression  is  third  order  in  the 
time  interval  At.   The  second  term  may  be  evaluated  by 
substitution  of  the  differential  equation  (2.3).   If 
(2.3)  is  differentiated  with  respect  to  the  independent 
variable  time,  then  -w,, ,  is  given  by 

(3.2)  =  .(|£  w,J,^  +  (Is  W,J,. 


-  [A(f,^+g,y)],^+  [B(f,^+  g'y^l'y  ' 

In  C3-2)  the  matrices  A  and  B  are  the  Jacobians  f.  ,  and  g. 
respectively  and  are  given  by 
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ACw)^ 


-1 


0 


^x2+  i^yS 


XY 


(7-3)X        (7-l)Y     (i-r) 


-Y 


-X 


0 


7XZ+(1-7)X(X^+Y^)  -7Z+  ^(3X^+Y^)   (7-l)XY      -jX 


(3. 5) 


B(w)^ 


0 


XY 


0 


-Y 


-1 


-X 


^  Y^  ^  X2 


0 


0 


(7-l)X 


(7-3)Y        (1-7) 


7YZ+(1-7)Y(X^+Y^)   (7-l)XY   -7Z+  ^^(  3Y^+X^)    -7Y 


The  abbreviations 

X  =  m/p  ,    Y  =  n/p  ,    Z  =  E/p 

have  been  used. 

Equations  (2.3)  and  (3.2)  when  substituted  into 
(3-1)  yield 


(5.3') 


w(t  +  At)  =  w(t)  +  At(f.   +  g.  )  + 


At' 


+  ^{[A(f,^+g,y)],^+[B(f,^+g,y)],y} 


To  preserve  the  high  order  accuracy  the  terms  containing 
first  derivatives  must  be  approximated  by  centered  differences 
while  the  term  containing  the  second  derivatives,,  since  it 
is  0    (At  ),  may  be  approximated  by  uncentered  differences. 
However,  to  center  the  entire  expression  consider  the  vectors 
C  evaluated  at  the  half  interval  (see  Fig.  la): 

C(AX/2)  =  fU^^y)+f(^^y)  +^(A++A)[f(x+,y)-f(x,y)  + 
(3.^a) 

+  -If[g(x^.y^)+g(x,y^)  -  g(x^,y~)-g(x,y")]]  . 

To  form  C(-AX/2),  we  wish  to  translate  the  argument  X  by 
the  amount  -Z\X.   Multiplication  by  the  translation  operator 
T  gives 

(X) 

{3.4b)  C(-AX/2)  =  T_^   C(AX/2)  . 

Similarly  in  the  y  direction 

C(AY/2)  =   gC^^y^)+g(x,y)  +  A^B++B)[^[f(x+,y+)+ 

(3-4c) 

+  f(x+,y)-f(x^,y~)-f(x,y")]+g(x,y+)-g(x,y)}, 

and 

(Y) 
(5.4d)  C(-AY/2)  =  T_^   C(AY/2)  . 

We  have  used  the  notation  A"^  =  A(X+AX:],  B"^  =  B(Y+AY),  etc. 
Then  the  difference  approximation  to  Eq.  (2.1)  is 

(•3.5)  w(t+At)  -  w(t)+A[[C(AX/2}-C(-AX/2)]+[C(AY/2)-C(-Ay/2)]}  . 
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Here  A  =  At/A  where  A  =  Ax  =  Ay.   It  Is  seen  then,  that  by 
centering  the  difference  scheme  has  one  of  the  properties 
of  the  differential  equation;  It  Is  In  divergence  free  form, 


-AX/2 


-AY/ 2 


AY/2 


AX/2 


Figure  la 

The  solution  of  the  axlally  symmetric  equations  can 
be  achieved  more  simply  by  a  two  step  calculation  analogous 
to  the  one  developed  by  R.  Rlchtmyer  [4].   Two  step  methods 
have  the  advantage  of  not  requiring  the  computation  of  the 
matrices  A  and  B  and  hence  matrix-vector  multiplication 
are  eliminated.   A  variant  of  the  two  step  Lax-Wendroff 
scheme  In  two  dimensions  Is: 


Step  1:   generates   temporary  data  at  the  four  points  (•) 
(k  =  1  +  1/2,  £   =    j   +   1/2,  t+At)  by  the  nine  points  (x) 
shown  In  Pig .  lb . 

w    (t+At)=  ^(w        +  w        ^+w   11+ 

k+-2,^+-2    k+  2'-^"  "2    ^~  2'-^+  '2 


(3.6) 


+  w 


k-  2,i-  2 


1^  ^1  ^^,  ,  1  „.  1  -f,_  1  .,  1  + 


k+  -^,£+  -|   k-  2,-^+  2 
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+  ^     1       1-^1       1 


)  +^[, 


=    1     1  'S    1  ^ 


+ 


+  g 


1  1  "  S   1 


l) 


6 — 

-^ 

• 

• 

• 

r 

• 

T 

— 3 

k- 

1^ 

^ 

Figure  lb 


Step  2: 


t+At 


(3-7) 


w.  .(t+At) 


=  w,_.Ct)  .^{f,^,_j-  f,_l,j+f     )^ 


Figure  lb ' 
The  notation  rj^^^^    expresses  the  fact  that  the  x-centered 

X 

difference   In   f  obtained   at   time   t+At   should  be   evaluated 
from  the    formula: 
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t+At 
^A    =    ^(^    1     1^  "  ^f^    1     1^  "^  ^^^.   1.1^" 

-  f(w_   ^  _   ^) 
and  similarly  for  g^   .   Fig.  (lb')  shows  as  solid  lines 

y 

difference  approximations  to  f,   and  g,   at  t+At  while  the 

•^      y 

dashed  lines  are  approximations  at  t.   The  one  dimensional 
version  of  this  scheme  was  originally  suggested  by 
B.  Wendroff.   The  method  centers  all  quantities  at  the 
point  (i,j,t+At/2)  by  the  averaging  procedure  contained 
in  Eq.  (3-7),  Just  as  is  done  by  the  method  of  Taylor 
expansion  in  Eq .  (5.I). 

For  axially-symmetric  flows  difference  equations  (3-6) 
and  (3-7)  are  modified  slightly.   By  inserting  the  terms 

k+-|,i+i  k+-i,^--i  k-^,£+^     k-  -|,i-  ^ 
on  the  right  hand  side  of  Eq .  (3.6)  and 

At  1  t+^t 

C3.7a)  |I  [i.^^.(t)+i^^^(t+At)],  i^^_g(t+At)=  ^(i_   ^  _   ^) 

_  "p '  J_  "2 
on  the  right  hand  side  of  Eq .  (3-7)5  the  two  step  method 
can  be  used  for  axially  symmetric  flows.   Equations  (3.6a) 
and  (3-7a.)  center  the  inhomogeneous  term  at  point  (i,j,  • 
t+  -p   At)  consistent  with  the  other  terms  in  Eqs .  (3-6)  and 
(3.7)- 

Equations  (3-6)  and  (3.7)  may  be  linearized  by  allowing 
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f  =  Aw  and  g  =  Bw  where  A  and  B  are  considered  constant. 

Then,  combining  equations  i^.6)    and  (3-7)  and  substituting 

Q  l(k  .eAx+k  mAy) 
w(iAx,mAy)  =  w  e         ^  ,    one  obtains  the  amplifica- 

tion matrix  [ 5] 

G'  =  I  +  iAA{-|sin  i    +\   sin^  cos  t\]    +  iBA[^  sin  ri  + 

2  2 
(3.8)     +  ^  Sinn  cos  ^  +  ^^    (cos  e-l)(l+cos  ri)  + 

+  ^^  (cos  n-l)(l+cos  U-  ^^t^^  A^sine  sinn  , 


where  ^  =  k  Ax,  tj  =  k  Ay  and  A  =  At/ Ax. 

It  may  be  verified  that  (3-8)  can  be  written  as 

G'    =   I+iA[A    sin   ^(1±£°^)+  B   sin   ti(5±^^)} 


(5-9) 


^    [A>y/[l-COS     4)  (1+COS    T])     +By(l-COS     Tl)Cl+COS     ^)] 


For  CjTI  ^"*  1  the  approximations 
Sin  1=1,    1-cos  i    =    -   4^/2  ,    1+cos  ^  =  2-^^/2  , 


may  be  used.   Equation  (3-9)  becomes 

2 
G'    =    I    +   iA(A4  +Bti)    -  ^    (A^  +Bii)2    +    <:^(e^,Tl^) 


Then,    modulo    terms    of   third    order   in   ^    and  r] , 

^,    ^    giA(A^+BTi) 

This  is  the  exact  amplification  matrix  of  the  differential 

equation  w,  +  Aw   +  Bw  .   Hence  it  is  clear  that  the  set  of 
t     X     y 

-      Ih      - 


equations  (3-6)  and  (3.7)  are  second  order  accurate. 

4 .   A  non-linear  Instability 

In  a  recent  paper  Lax  and  Wendroff  presented  a 
stability  proof  of  the  second  order  difference  method 
given  by  Eq .  (5-5) •   We  call  the  amplification  matrix 
associated  with  the  linearized  form  of  Eq .  (5-5)  G; 
it  is  given  by 


2 
Gii,r\)    =  I  +  1  (A  sin  ^  +  B  sin  t^  )  -  A  (l-cos  |] 


(4.1) 


T.2r-,       ,    ,AB+BAv   .   . 

B  (l-cos  ri)  -  ( 5 }  sm  t,    sm  ri  . 


In  order  to  prove  stability  one  has  to  show  that  all 
eigenvalues  of  the  amplification  matrix  are  ^  1  in 
absolute  value  (von  Neumann  condition);  this  follows 
if  the  inner  product 

(^•1)'  (Gq  ,  q)  <  1 

for  all  unit  vectors  q.   In  fact  (4.1)'  implies  that 

I  n  I 

|G  I  _<  K  for  Integer  powers  of  n,  i.e.,  n  =  1,2,...  . 

The  difference  scheme  associated  with  G  will  be  stable, 

at  least  if  the  coefficients  of  G  are  constant,   since 


* 


Recently,  Lax  and  Nirenberg  have  shown  how  to  extend 
this  conclusion  to  G  with  variable  coefficients. 
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the  solutions  of  the  difference  scheme  (obtained  by  n 
successive  multiplications  of  the  difference  operator) 
will  be  bounded.   This  is  a  consequence  of  the  isometry 
of  the  Lp  norm  when  a  fourier  transform  is  made  on  the 
linear  difference  operator  to  obtain  the  amplification 
matrix.   An  analysis  presented  in  [2]  gives  the 
following  inequality: 

KGq,  q)|^  <  1  -  Uql^Cl-  8 | Aq | ^] (l-cos  ^)^  - 


(4.2) 


-  |Bq|^(l-  8|Bq|^)(l-cos  x\)^ 


This  expression  is  less  than  unity  if  the  requirements 

A^  <  1/8  ,       B^  <  1/8  , 

are  satisfied.   For  the  case  of  equal  spatial  stepsize 
Ax  =  Ay  =  A   this  condition  can  be  expressed  as 

(4.3)   ^—  <  —zi ,  a     =   max  eigvalue  of  A  or  B. 

^    a  /8 

In  [1]  is  was  shown  empirically  that  the  above  stability 
condition  is  too  stringent;  it  can  be  exceeded  in  practical 
applications  without  causing  instability;  in  particular 
At/A  =  l/j/2a,  which  is  twice  as  large  as  permitted  by  (4.3) 
yielded  stable  results.   It  was  even  possible  to  go 
slightly  (^  lO/o).  beyond  this  numerical  estimate.   However, 
for  the  class  of  problems  considered  in  this  paper,  even 
satisfying  l/lO  of" condition  (4.3)  led  to  numerical 
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Instability  during  the  early  transient  portion  of  the 
flow  (pointwise  values  of  the  density  p  and  energy  E. 
>  10-^  ~  10^  after  which  the  density  became  negative)  . 
During  this  time,  the  shock,  having  been  produced 
impulsively  at  the  face  of  the  body,  was  in  the  process 
of  detaching  itself;  instabilities  were  observed  in 
two  regions:   within  the  shock  close  to  the  axis  of 
symmetry  and  at  the  corner  of  the  body.   The  corner 
problem  could  be  removed  if  a  first  order  difference 
approximation  (discussed  in  Section  7)  was  used  on  the 
top  of  the  body.   But  within  the  shock  negative 
densities  were  still  obtained. 


supersonic 
flow 


supersonic 
flow 


body 


Figure  Ic 

Figure  2  shows  schematically  the  growth  of  the 
solution,  at  the  shock,  when  an  instability  is  present. 
The  highly  oscillatory  nature  of  the  solution  will  grow 
in  time  if  the  calculation  is  allowed  to  continue.   The 
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data  given  in  Figures  2  and  3  were  obtained  from  results 
which  were  an  attempt  to  eliminate  the  early  instability 
(after  approximately  6o  cycles)  that  was  obtained  with 
Eq.  (3.5).   Here,  a  modification  in  the  coefficient  of 
the  second  term  [5]  in  {^ Aa) ,    A/4(A^+A) -D  and  (5-^c), 
A/4(b'^+B)-D  where  D  >  1  allowed  the  calculation  to  be 
carried  out  as  far  as  5OO  cycles  but  eventually  the 
phenomenon  shown  in  Figure  2  resulted.   For  D  >  1,  the 
accuracy  of  the  difference  scheme  is  no  longer  second 
order.   The  result  of  a  case  for  D  =  3  is  given  in 
Figure  3  which  shows  the  shock  position  for  successive 
times.   At  the  point  where  the  shock  position  becomes 
stationary  (at  approximately  490  cycles),  however, 
instability  results.   The  growth  of  the  amplitude  of 
short  wave  length  components  of  the  solution  for  density 
p,  is  shown  in  Figure  2.   This  was  localized  near  the 
axis  of  symmetry  at  the  shock. 

The  fourth  order  artificial  viscosity  term  which  is 
discussed  in  [2]  was  also  tested  with  Eq.  (3.5)-   It  had 
virtually  no  effect  on  the  stability  of  the  numerical 
results . 

These  results  were  surprising  but  can  be  explained 
by  the  following  argument.   Inequality  (4.2)  together  with 
the  stability  condition  show  that  the  eigenvalues  of  G  are, 
for  ^,7]   ^  0,    definitely  less  than  one;  this  means  that  high 
amplitude  waves  are  damped.   However,  if  |Aq|  or  |  Bq  |  happen 


to  be  zero,  this  damping  is  absent.   Now  |Aq|,  |Bq|  can 
vanish  if  and  only  if  zero  is  an  eigenvalue  of  A  or  B . 
The  eigenvalues  of  A,  B  are  u  +  c,u  and  v+c,v.   Hence  for 
points,  lines  or  surfaces  in  the  flow  where 

Condition  I:   u/c  =  1,  or  v/c  =  1, 

(c  is  the  local  sound  speed)  or 

Condition  II:   u  =  0  =  v, 

one  of  the  eigenvalues  becomes  zero.   Condition  I 
corresponds  to  local  one  dimensional  sonic  flow  while 
condition  II  corresponds  to  flow  at  stagnation  points. 
If  the  sonic  line  is  extended,  it  will  pass  within  the 
shock  (see  Fig.  Ic)  so  that  near  the  axis  of  symmetry 
since  the  flow  is  locally  one-dimensional,  t]  =  0, 
Condition  I  occurs;  i.e.,  there  will  be  no  damping  for 
the  associated  linear  difference  scheme  (4.2)  and  hence 
I (Gq  ,  q) I   =1.   Since  this  condition  of  neutral  stability 
occurs  in  a  part  of  the  flow  where  the  basis  of  applying 
linear  analyses  are  weakest,  (where  flow  properties  vary 
most  rapidly)  a  nonlinear  statement  of  stability  is 
suspect.   Such  second  order  schemes  (including  the  two 
step  method  given  by  Eqs .  (3.6)  and  (3-7))  will  probably 
need  additional  damping  terms  for  successful  results. 
These  terms  must  be  third  order  so  as  to  preserve  the 
numerical  accuracy. 
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-,     I  a.  (co)  -  a.  (co'  )  I 

(5.3a)  a^  =  -i  '^   _   ^ — : —  ,  1=  1,2,3,  j  7^1,  k^l,J. 

(a.  -a  . )  (a.  -a,  ) 

Here  the  barred  quantities  are  understood  to  be  averages 
over  the  Interval  (a),co'),  and  «■ .    may  be  a  function  of  a 
or  constant.   This  will  insure  that  if  the  flow  is  locally 
constant  the  a.  will  vanish.   A  similar  relationship  can 
be  written  for  p.  in  terms  of  a(B) .   With  this  choice  for 
a   and  ^,    the  coefficients  of  Eq.  (5.2)  can  be  put  in  the 
form 

1=1, j^i   ^ 

^  :^ 

(5- 3b)  a  =  -  2Z  ct.(a  +  a  )  , 

1    i=l7j^i  ^  J   ^ 

3 


ap  =  ^  a. 

^    i=l   ^ 


Then  the  matrices  defined  in  (5-2)  can  be  given  by 

Q"^  =  ZI   a,[(A-a,I)(A-a,  I)  , 
k7^1,j 
(5.4) 

3 


Q^  =  21   p.(B-aI)(B-a^I)} 
1=1,  jVi 
k?^i,  J 


It  is  clear  that  with  such  a  choice  of  coefficients  the 
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eigenvalues  of  Q^  or  Q'^  will  be  proportional  to  the 
difference  of  the  eigenvalues  of  A  or  B  in  the  interval 
(cd,cd'),  i.e. 

Hence  since  there  will  always  be  at  least  one  non-zero 
eigenvalue,  the  viscosity  will  contribute  to  the  solution 
in  regions  where  the  solution  varies  most  rapidly.   This 
will  insure  us  that  the  difference  scheme  will  have 
dissipation. 

Hence  the  complete  difference  scheme  obtained  by 
combining  Eq .  C3.5)  with  the  artificial  viscosity  given 
by  Eq.  (5.I)  is 


w(t+At)  =  w(t)  +  A{D  [L  C(x)  +^  Q^'w]  + 


(5.5) 


+  D  [L  C(y)  +  -^  Q'^d'w]} 


We  have  used  L  C(x)  =  -^  { C(x+Ax) +CCx) }  =  CCAx/2)  as  the 

X  c- 

I 

forward  averaging  operator,  D  and  D  as  the  centered  and 

X,         X 

forward  difference  operators  respectively. 

For  the  Lagrangian  system  of  hydrodynamic  equations 
these  coefficients  reduce  to  that  given  by  Lax  and 
Wendroff  i.e. , 


^0 

- 

0 

= 

^1    ' 

a^ 

= 

1 
2 

K 

C    CO 

-c(a) ' 

) 

52 
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In  the  next  section  we  use  this  Lagranglan  form  of 
the  viscosity  to  show  stability  of  the  Lax-Wendroff 
scheme  with  the  added  artificial  viscosity. 

6 .   Stability  analysis  of  second  order  scheme  with  viscosity 

In  order  to  study  the  stability  of  the  difference 
schemes  we  first  linearize  the  equations  by  taking  the 
matrices  A  and  B  to  be  locally  constant.   We  assume  that 
the  real  matrices  A,  B  are  symmetric  or  at  least  they  can 
be  symmetrized  by  a  similarity  transformation.   Also,  for 
simplicity,  let  Z\x  =  Ay  =  At  =  1 .   We  also  take  the 
viscosity  to  be  of  the  form 

(6.1)      Q(A,B)  =  Q^'CA)  +  Q^(B)  =  |  a  A^  +  1  p  B^  , 

where  a  and  p  are  constant.   This  choice  of  Q  is  then 


substituted  into  (5.I)  and  the  result  is  linearized, 

y 


1  22  1        22 

l.e.T;aA     D     w+ttRB     D     w,    where   the   second   difference 

2  X     2  "^      y 


2 
operator  is  represented  by  D  .   Take  the  Fourier  transform 

of  the  above  and  add  it  to  the  amplification  matrls  given 

by  Eq .  (4.1);  the  result  is: 

Gii,r])    =  I  +  i(A  sin  ^  +  B  sin  t^)  -  (1  +  -|  a)A^(l-cos  ^) 

(6.2) 

,-,  ,  1  „  ^-o2,-,         ,    AB  +  BA   .   ,. 
-  (1  +  -p  Pi)B  (1  -  cos  T] )  -  p sm  ^  sm  -q  . 

The  proof  basically  parallels  that  presented  by  Lax  and 
Wendroff  [2];  we  give  all  steps  for  completeness. 
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We  wish  to  estimate  the  quantity  (Gq  ,  q)  for  all 
unit  vectors  q.   G  is  separated  into  its  real  and 
imaginary  parts 

G  -  X  +  iY 
where 

Y  =  Asln4+BsinTi 
and 

X  =  I  -  R 
with 

(6.3)   R  =  (1+  "I  a)A^&    +  Cl+  ^  P)B^  +  M^  sin  ^  sin  t^ 

We  have  used  the  definitions  (h)  =  1  -  cos  i   and 
5=1-  cos  r\ .      Now 

CGq  ,q)  =  (Xq  ,  q)  +  l(Yq  ,  q)  . 

Since  A  and  B  are  real  and  symmetric,  it  follows  that  X 
and  Y  are  real  and  symmetric;  hence 

C6.4)         |(Gq,  q)|^  =  (Xq,  q)^  +  (Yq  ,  q)^  . 

R  can  be  given  by  the  identity 

(6.5)   R  =  I  a2@2  +  1  B^2  +  1  y2  +  1  ^  a2(5)  +  I  pB^  . 

The  proof  proceeds  by  computing  the  right  hand  side  of 
Eq.  (6.4).   Using  the  abbreviations  x  =    (Xq  ,  q), 
y  =  (Yq  ,  q),  and  r  =  (Rq  ,  q)  the  absolute  value  of  the 
real  part  of  the  inner  product  of  the  amplification 
matrix  is 
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(6.6)  x^  =  (1  -  r)^  =  1  -  2r  -  r^  . 

Now  the  Identity  In  R,  Eq .  (6.5),  Is  used  to  obtain  the 
second  term  in  Eq.  (6.6)    while  the  definition  of  R,  given 
by  Eq .  (6.3) >  is  used  for  the  evaluation  of  r  .   Taking 

the  Inner  product  of  Eq .  (6-5)  and  using  the  Schwartz 

2        I    I  2     2 
Inequality  (Yq  ,q)  =  |Yq|   ^  y  ,  we  obtain 

(6.7)  r  <  -|  a2(h)2+  1  B^2  +  1  y2  _^  1  ^^2^  ^  1  ^^2-^  _ 

2 
To  obtain  | (Rq  ,q)|  ,  multiply  Eq .  (6.5)  by  q  and  form 

the  Inner  product.   One  of  the  resulting  terms, 

Re(Aq  ,  Bq)  sin  ^  sin  t]  ,  may  be  estimated  by  the  Schwartz 

inequality 

j<  I  Aq|  I  sin  ^  I  |Bq|  |  sin  r]  |  _<  [  |  Aq  |  (h)  +  |Bq  |  "%] 

1     2 
The  trigonometric  inequalities  (S)  >  ^   sin  i   and 

J  _>  -p  sin  ^   have  been  used.   Using  this  result,  Eq.  (6.7) 

becomes 

r<  iAq|^  (2  +^a)(S)+  |Bq|^  (2  +  ^  ^)J 


or 


r   < 
(6.8) 


2  -  (|Aq|2)^   ^[(2  +  ^)2  +  (2+|)(2+|)} 


+  ClBql^)^  [(2+|)2  +  (2+|)(2+|)}  . 
Substituting  Eqs .  (6.8)  and  (6.7)  into  Eq.  (6.6)    gives 
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x^  +  y^  <  1  -  |Aq|^{{l  +  a)  -  a^|Aq|^}@^ 

C6.9) 


-  |Bq|2[(i+^)-pllBq|23$2  . 


We  have  used  the  trigonometric  inequalities  (for  4,'n'v7r) 
(h)   _>  (h)  and  ^  _>  ^  (Short  wave  length  disturbances  must 
be  damped  so  that  their  growth  will  not  affect  the  numeri- 
cal solution;  see  Figure  2)   to  obtain  Eq .  (6.9).   The 
right  hand  side  is  less  than  one  if 

(6.10)   lAql^  =  |a(A)q|2^i^  ;   |Bq|2=  |a(B)q|2^^  . 

a   and  p   are  abbreviations  for  the  quantities  in  brackets 
in  Eq.  (6.8) . 

It  is  interesting  to  note  that  Eq .  (6.9)  is  of  ■  .;■ 
the  same  form  as  Eq.  (4.2).   The  argument  presented  in 
Section  4  on  nonlinear  stability  would  appear  to  hold 
in  the  present  case.   It  must  be  concluded  (in  light  of 
the  numerical  experiments  performed  with  Eq .  (5.5))  that 
linearization  eliminates  many  important  features  of  the 
difference  scheme,  e.g.  nonlinear  terms  (see  Richtmyer 
and  Morton  [7])-   For  the  simple  case  of  equal  coefficients 
for  the  viscosity,  Eq .  (6.10)'  becomes 

.;       a2(A)  <        ^  +  g    =  f(a)  . 

2(2  +  |]2 

The  max  of  f(a}.  is  for  a  =  2,  f(2]  =  1/6.   Hence  there  is 
hope  that  the  scheme  will  be  stable  for 
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C6.li)  ii"t*'i-^   ' 

while  a  similar  expression  holds  for  a(B) .   This  is 
somewhat  larger  than  the  value  1/|/B  (a=  0),  Lax  and 
Wendroff  obtained.    The  term  given  by  Eqs .  (5.I)  and 
(6.1)  which  is  an  artificial  viscosity,  clearly  has  a 
stabilizing  influence  on  the  difference  scheme  Eq .  (5-5) > 
since  the  allowable  time  increment  is  larger.   Similar 
results  hold  for  the  more  complicated  forms  of  viscosity 
used  in  the  numerical  calculation,  i.e.  Eqs.  (5.I)  and 
(5-2). 

The  test  of  the  stability  condition  was  made  using 
a  constant  to  represent  the  right  hand  side  of  (6.11). 
This  numerical  value  was  fed  into  the  machine  via  a  data 
card  and  the  resulting  solution  was  checked  every  30 
cycles  of  computation.   No  other  changes  were  made  and 
the  results  appear  in  Table  1.   The  theoretical  limit 
r  <  .408  (which  is  close  to  the  Lax  Wendroff  value  of 
•55^)  agrees  within  the  accuracy  obtained  with  the 
numerical  estimate  of  Table  1.   Additional  experiments 
will  be  performed,  especially  for  the  case  of  axially 
symmetric  flows. 


Although  a  and  p  are  not  arbitrary  constants  but  depend 
on  the  solution  about  which  linearization  is  performed, 
locally  the  stability  limit  imposed  by  Eq.  (4.3)  may  be 
relaxed  for  the  condition  given  by  Eq.  (6.11).   In 
practical  applications^,  however,  the  difference  is 
unimportant  and  Eq.  (4.3)  may  be  used. 
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Table  1 


r  = 

a 
max 

At 
A 

Condition 

.65 

unstable  at  90  cycles 
(negative  densities) 

•55 

negative  velocities  at 
cycles  near  stagnation 

350 

point 

.45 

negative  velocities  at 
cycles  near  stagnation 

640 
point 

•35 

no  negative  densities  ( 
velocities  at  25OO  eye 

Dr 

Les  . 

The  fact  that  the  velocities  became  negative  for  some  values 
of  r,    is  in  itself  not  an  instability.   The  solution  in  the 
region  of  the  stagnation  point  did  start  to  diverge  from 
the  known  solution  when  negative  velocities  were  observed. 
As  the  computation  proceeded  the  situation  became  worse 
since  the  negative  values  of  the  x  component  of  velocity 
propagated  into  the  mesh  away  from  the  stagnation  point. 
At  this  point  the  calculation  for  the  case  was  terminated. 
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7 .   Numerical  results 

The  results  of  two  calculations  are  presented.   The 
difference  scheme  (5.  5),  which  indudes  Hie  artificial  viscosity 
(5.1),  was  used  to  obtain  these  results.   The  coefficient 
IC    in  Eq.  (5.3a)  was  initially  set  equal  to  2  but  during 
the  course  of  the  calculation  it  was  reduced  to  1.   The 
free  stream  Mach  number  is  4.3  and  j,    the  ratio  of 
specific  heats  is  1.4.   The  flow  was  started  Impulsively 
by  using  the  free  stream  conditions  as  initial  data.   In 
Figure  4  we  plot  stagnation  density  vs.  number  of  cycles. 
Beyond  the  2,340-th  cycle  the  values  of  density  everywhere 
in  the  flow  field  are  constant  to  four  figures;  clearly 
the  method  converges.   The  accuracy  in  the  calculation  of 
the  stagnation  density  is  exceptionally  good  while  the 
error  of  .96/0  in  the  computed  stagnation  pressure  is 
consistent  for  the  method.   Figure  5  shows  the  density 
for  values  of  y  as  a  function  of  x.   The  distance  y  has 
been  normalized  by  the  radius  of  the  body.   The  corner  of 
the  body  is  at  x  =  1.28.   For  y  =  1,  it  is  seen  that  there 
is  a  rapid  variation  of  the  density  corresponding  to  the 
corner  expansion.   For  y  =  2.2,  the  density,  after  the 
shock,  decreases  linearly  corresponding  to  the  fluid 
expanding  past  the  body.   The  complete  flow  field  is 
shown  in  Figure  6.   Here  the  position  of  the  shock  and 
sonic  line  are  shown  with  respect  to  the  body.   The 
dashed  curve  represents  the  sonic  line  computed  with 
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the  scheme  of  Lax  given  in  1952. 

2      P 

(7.1)      w(t+At)  =  w(t)  +  ACD  w  +  D  w)  +  (D  w  +  D  w)  . 

X    y      ^    y 

Equation  C7-1)  is  a  first  order  method  and  the  point  of 
attachment  of  the  sonic  line  to  the  body  is  greater  than 
3  space  increments  from  the  body  corner.   On  the  other 
hand  the  high  order  scheme,  Eq.  (5.5),  gives  the  position 
of  the  sonic  line  correctly  i.e.  within  a  fraction  of  one 
mesh  point  of  the  corner. 

As  shown  in  Figure  7?  the  Mach  number  at  the  corner 
is  0.975.   One  mesh  point  to  the  right  on  top  of  the  body 
the  Mach  number  is  2.04  while  it  is  0.666  one  mesh  point 
below  the  corner,  on  the  body  face.   Beyond  the  corner 
the  Mach  number  increases  very  rapidly,  indicating  a  strong 
rarefaction  wave. 

Thus  the  difference  scheme  is  able  to  handle  such 
rapid  variations  of  the  flow  quantities;  this  is  again 
indicated  by  the  distribution  of  pressure  ratio  along 
the  body  given  in  Figure  8.   In  Figure  9  the  entropy 
distribution  along  the  stagnation  streamline  is  shown. 
The  theoretical  value  of  entropy  at  the  stagnation  point 
can  be  calculated  exactly;  it  is  2.24  while  the  computed 
value  is  2.20.   The  abrupt  entropy  change  near  the 
corner  and  the  oscillation  on  top  of  the  body  are  very 
likely  due  to  the  fact  that  the  vertical  component  of 
velocity  was  computed  on  the  top  of  the  body  rather  than 
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prescribed  to  be  zero.   This  was  done  so  as  to  observe  the 
number  of  space  steps  needed  to  approach  the  correct  value. 
This  meant  that  improper  values  of  entropy,  mass,  etc.,  were 
introduced  into  the  flow  because  of  false  boundary  conditions. 
The  entropy  distribution  A  above  the  body  is  shown  and  it 
approaches  the  theoretical  correct  value.  The  sharp  rise  in 
the  value  of  the  entropy  at  the  corner  is  presumably  due  to 
the  artificial  viscosity  acting  on  the  steep  rarefaction 
gradients  at  the  corner.   In  order  to  resolve  this  problem, 
additional  work  is  necessary.   Figures  10  and  11  show  some 
results  for  the  case  where  M  =^  10  and  7  =  1.17-   In  this 
calculation  it  was  found  necessary  to  introduce  the  function 

ic  =  K^[i+exp  (i-ia^iA^)}      if  |a^|A^_<i, 

or 

'^  =  1  otherwise; 

for  the  coefficient  of  viscosity  given  in  Eqs .  C5-3a).   '^-i 
was  chosen  to  be  of  the  order  of  the  velocity  behind  the 
strong  normal  shock,  while  '<^„  was  between  one  half  and  one; 
|u  I  is  the  magnitude  of  the  velocity.   This  has  the  effect 
of  providing  additional  damping  near  the  stagnation  point. 
This  damping  was  not  required  for  the  previous  case  (M  =  4.5) 
since  the  shock  was  well  away  from  the  body.   It  is  felt,  for 
the  case  M  =  10,  that  if  unequal  spacing  were  chosen  in  the 
X  direction,  so  that  there  would  be  more  steps  between  the 
shock  and  the  body,  additional  damping  would  be  unnecessary. 
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350 
2  80 

210 
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70 


SONIC  LINE  N=420  CYCLES 


Time  Dependent  Solution  Showing  Shock  Position  for  Modified   Lax-Wendroff 
Scheme.    Numbers  Refer  to  Computational    Cycles. 


Figure  3 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Commission,  nor  any 
person  acting  on  behalf  of  the  Commission: 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information  contained  in  this  report, 
or  that  the  use  of  any  information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus, 
method,  or  process  disclosed  in  this 
report . 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Commission"  includes  any  employee  or 
contractor  of  the  Commission,  or  employee  of 
such  contractor,  to  the  extent  that  such  em- 
ployee or  contractor  of  the  Commission,  or 
employee  of  such  contractor  prepares,  dis- 
seminates, or  provides  access  to,  any  infor- 
mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with 
such  contractor. 


